home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1997 August / Walnut Creek CDROM.7z / VOL_400 / 446_01 / DOC / DPSTART / TEST2 / WAVE.C < prev   
Encoding:
C/C++ Source or Header  |  1996-04-18  |  3.2 KB  |  87 lines
  1. #include <Arrays_real.h>
  2. #include <CurvPlot.h>
  3.  
  4. // forward declarations:
  5. void scheme (Vec(real)& up1, Vec(real)& u, Vec(real)& um1, real tstop, real C);
  6. void setInitialCondition (Vec(real)& u0, Vec(real)& um1, real C);
  7. void plotSolution (Vec(real)& u, CurvPlotFile& plotfile, real t, real C);
  8.  
  9. int main (int nargs, const char** args)
  10. {
  11.   initDIFFPACK (nargs, args);
  12.   cout << "Give number of intervals in (0,1): ";
  13.   int i; cin >> i;
  14.   int n = i+1; // number of points;
  15.   Vec(real) up1 (n);  // u at time level k+1
  16.   Vec(real) u   (n);  // u at time level k
  17.   Vec(real) um1 (n);  // u at time level k-1
  18.   cout << "Give Courant number: ";
  19.   real C; cin >> C;
  20.   cout << "Compute u(x,t) for t <= tstop, where tstop = ";
  21.   real tstop; cin >> tstop;
  22.  
  23.   setInitialCondition (u, um1, C);  // initialization
  24.   scheme (up1, u, um1, tstop, C);   // finite difference scheme
  25.   return 0;
  26. }
  27.  
  28. void scheme (Vec(real)& up1, Vec(real)& u, Vec(real)& um1, real tstop, real C)
  29. {
  30.   int  n = u.size();   // length of the vector u (no of grid points)
  31.   real h = 1.0/(n-1);  // length of grid intervals
  32.   real dt = C*h;       // time step, assumes UNIT wave VELOCITY!
  33.   real t = 0;          // time
  34.   CurvPlotFile plotfile(CaseName); // "databank" (file) with all the plots
  35.  
  36.   plotSolution (u, plotfile, t, C); // plot the initial string displacement
  37.   int  i;                           // loop counter
  38.   int  step_no = 0;                 // current step number
  39.  
  40.   while (t <= tstop)
  41.     {
  42.       t += dt;  step_no++;
  43.       for (i = 2; i <= n-1; i++)
  44.         up1(i) = 2*u(i) - um1(i) + sqr(C) * (u(i+1) - 2*u(i) + u(i-1));
  45.       up1(1) = 0;  up1(n) = 0; // boundary condition at x=0 and x=1
  46.       plotSolution (up1, plotfile, t, C);
  47.       um1 = u;  u = up1;       // update data structures for next step
  48.       if (step_no % 100 == 0)  // write a message for every 100th step:
  49.         cout << oform("time step %4d: u(x,t=%6.3f) is computed.\n",step_no,t)
  50.              << flush; // the flush forces immediate output
  51.     }
  52. }
  53.  
  54. void setInitialCondition (Vec(real)& u0, Vec(real)& um1, real C)
  55. {
  56.   int  n = u0.size();  // length of the vector u
  57.   real x;              // coordinate of a grid point
  58.   real h = 1.0/(n-1);  // length of grid intervals
  59.   real umax = 0.05;    // max string displacement
  60.  
  61.   for (int i = 1; i <= n; i++)  // set the initial displacement of the string
  62.     {
  63.       x = (i-1)*h;
  64.       if (x < 0.7)
  65.         u0(i) = (umax/0.7) * x;
  66.       else
  67.         u0(i) = (umax/0.3) * (1 - x);
  68.     }
  69.   for (i = 2; i <= n-1; i++) // set the help variable um1:
  70.     um1(i) = u0(i) + 0.5*sqr(C) * (u0(i+1) - 2*u0(i) + u0(i-1));
  71.   um1(1) = 0; um1(n) = 0;  // dummy values, will not be used in the scheme
  72. }
  73.  
  74. void plotSolution (Vec(real)& u, CurvPlotFile& plotfile, real t, real C)
  75. {
  76.   int  n = u.size();       // the number of unknowns
  77.   real h = 1.0/(n-1);      // length of grid intervals
  78.   CurvPlot plot (plotfile);
  79.   plot.initPair ("displacement",               // plot title
  80.                  oform("u(x,%g)",t),           // name of function
  81.                  "x",                          // name of independent variable
  82.                  oform("C=%g, h=%g",C,h));     // comment
  83.   for (int i = 1; i <= n; i++)
  84.     plot.addPair (h*(i-1) /* x-value */, u(i) /* function value */);
  85.   plot.finish();
  86. }
  87.